function [Q,F,S,empty] = drawQ_ARW(restr,irs,phi,opt,W)
% Function draws Q from the space of orthonormal matrices satisfying the
% equality and sign restrictions using the algorithm in Arias,
% Rubio-Ramirez and Waggoner (2018) (modified to impose sign normalisation 
% in GK18).
% Inputs are:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run (cumulative) IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options
% - W: cell array of standard normal random variables. Use the following
% code to draw this before calling drawQ_ARW.m:

% W=cell(N,1);
% for jj=1:N
%     W{jj}=randn(N-(jj-1+f(jj)),N);
% end

eqRestr = restr.eqRestr;
signRestr = restr.signRestr;
f = restr.f;
s = restr.s;
vma = irs.vmaGK;
lrcir = irs.lrcirGK;
Sigmatr = phi.Sigmatr;
Sigmatrinv = phi.Sigmatrinv;
B = phi.B;
p = opt.p;
const = opt.const;
signNorm = opt.signNorm;                         
                            
N = size(Sigmatr,1); % Number of variables in VAR
mZ = size(eqRestr,1); % Number of equality restrictions
mS = size(signRestr,1); % Number of sign restrictions

% Construct matrix representing equality restrictions.
% F(phi) = [F_1(phi); ...; F_N(phi)], where F_i(phi) represent equality 
% restrictions on the ith column of Q, so F_i(phi)*q_i = 0.

if mZ == 0 % If no equality restrictions

    F = [];

else % If equality restrictions

    F = zeros(mZ,N);

    for ii = 1:mZ

        if eqRestr(ii,4) == 1 % Restriction on A0

            F(ii,:) = Sigmatrinv(:,eqRestr(ii,2))';

        elseif eqRestr(ii,4) == 2 % Restriction on A0^(-1)

            F(ii,:) = Sigmatr(eqRestr(ii,1),:);

        elseif eqRestr(ii,4) == 3 % Restriction on A_l
            
            B = reshape(B,[N*p+const,N]);
            Bl = B(const+(eqRestr(ii,3)-1)*N+1:const+(eqRestr(ii,3))*N,...
                eqRestr(ii,2));
            F(ii,:) = (Sigmatrinv*Bl)';
            
        elseif eqRestr(ii,4) == 4 % Restriction on LRCIR

            F(ii,:) = lrcir(eqRestr(ii,1),:);

        end

    end

end

% Order rows of F in terms of the column of Q restricted.

Ftilde = [];

if isempty(F)
    
    Ftilde = F;
    
else

    for jj = 1:N % For each column of Q

       % Find rows of F restricting jth column of Q.
       Fj = F(eqRestr(:,1)==jj & (eqRestr(:,4)==1 | eqRestr(:,4)==3)...
           | eqRestr(:,2)==jj & (eqRestr(:,4)==2 | eqRestr(:,4)==4),:);
       Ftilde = [Ftilde; Fj];

    end
    
end

% Construct matrix representing inequality restrictions (excluding 
% normalisations).
% S(phi) = [S_1(phi); ...; S_N(phi)], where S_i(phi) represent inequality 
% restrictions on the ith column of Q, so S_i(phi)*q_i >= 0.

S = zeros(mS,N);

for ii = 1:mS
    
    if signRestr(ii,5) == 1 % Sign restriction on IRF
        
        S(ii,:) = vma(signRestr(ii,1),:,signRestr(ii,3)+1)*signRestr(ii,4);
        
    elseif signRestr(ii,5) == 2 % Sign restriction on A0
        
        S(ii,:) = Sigmatrinv(:,signRestr(ii,2))'*signRestr(ii,4);
        
    end
    
end

% Order rows of S in terms of the column of Q restricted.

Stilde = [];

if isempty(S) 
    
    Stilde = S;
    
else
    
    for jj = 1:N % For each column of Q

       % Find rows of S restricting jth column of Q.
       Sj = S((signRestr(:,1)==jj & signRestr(:,5)==2) | ...
           (signRestr(:,2)==jj & signRestr(:,5)==1),:);
       Stilde = [Stilde; Sj];

    end
    
end

iter = 0;
empty = 1;

while iter <= maxDraw && empty == 1

    iter = iter + 1;

    qtilde = zeros(N);

    s1 = N - f(1);
    x1 = randn(N-f(1),1);
    w1 = x1./norm(x1);
    M1_tilde = [F(1:f(1),:); W{1}];
    [K,~] = qr(M1_tilde');
    K1 = K(:,N-s1+1:N);
    qtilde(:,1) = K1*w1;

    for jj = 2:N

        sj = N+1-jj-f(jj);
        xj = randn(sj,1);
        wj = xj./norm(xj);
        Mj_tilde = [qtilde(:,1:jj-1)'; ...
            F(sum(f(1:jj-1))+1:sum(f(1:jj)),:); W{jj}];
        [K,~] = qr(Mj_tilde');
        Kj = K(:,N-sj+1:N);
        qtilde(:,jj) = Kj*wj;

    end

    Q = zeros(N);

    for jj = 1:N

        if signNorm == 0

            % Normalise diagonal elements of A0 to be positive.

            if Sigmatrinv(:,jj)'*qtilde(:,jj) == 0

                % Randomly set sign of q_j to -1 or 1.
                Q(:,jj) = (1-2*(rand > 0.5))*...
                    (qtilde(:,jj)./norm(qtilde(:,jj)));

            else

                Q(:,jj) = sign(Sigmatrinv(:,jj)'*qtilde(:,jj))*...
                    (qtilde(:,jj)./norm(qtilde(:,jj)));

            end 

        else % Normalise diagonal elements of A0^(-1) to be positive

           if Sigmatr(jj,:)*qtilde(:,jj) == 0

                % Randomly set sign of q_j to -1 or 1.
                Q(:,jj) = (1-2*(rand > 0.5))*...
                    (qtilde(:,jj)./norm(qtilde(:,jj)));

           else

               Q(:,jj) = sign(Sigmatr(jj,:)*qtilde(:,jj))*...
                   (qtilde(:,jj)./norm(qtilde(:,jj)));

           end

        end

    end
    
    % Check whether proposed draw satisfies sign restrictions (if any).

    if mS > 0  % If sign restrictions

        signCheck = zeros(mS,1);
        restrCount = 0;

        for ii = 1:N
            
            for jj = 1:s(ii)
                
                restrCount = restrCount + 1;

                signCheck(restrCount) = S(restrCount,:)*Q(:,ii);
                
            end
            
        end

        % If IRFs satisfy sign restrictions, terminate while loop and 
        % record that set is non-empty.
        empty = (min(signCheck) < 0);
        
    else % If no sign restrictions
        
        empty = 0;
        
    end
    
end

F = Ftilde; % Replace F with reordered F
S = Stilde; % Replace S with reordered S

end